Monte Carlo simulation of granular fluids 
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Abstract 



An overview of recent work on Monte Carlo simulations of a granular binary mixture is presented. The results 
are obtained numerically solving the Enskog equation for inelastic hard-spheres by means of an extension of the well- 
known direct Monte Carlo simulation (DSMC) method. The homogeneous cooling state and the stationary state 
reached using the Gaussian thermostat are considered. The temperature ratio, the fourth velocity moments and 
the velocity distribution functions are obtained for both cases. The shear viscosity characterizing the momentum 
transport in the thermostatted case is calculated as well. The simulation results are compared with analytical 
predictions showing an excellent agreement. 
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c/> !1 Introduction 

]An usual way of capturing the dissipative nature of granular media is through an idealized fluid of smooth, inelastic 
L^^'hard spheres. Despite the simplicity of the model, it has been shown to be quite useful in describing the dynamics of 

Agranular materials under rapid flow conditions ^ |2] ■ The essential difference from ordinary fluids is the absence of 

energy conservation, leading to both obvious and subtle modifications of the Navier-Stokes hydrodynamic equations. 

In the context of kinetic theory, the Boltzmann and Enskog equations have been conveniently modified to account for 
I inelastic binary collisions. These kinetic equations have been used to derive the corresponding fluid dynamic equations 
i— H with explicit expressions for the transport coefficients for a monocomponent system |S] , and for a binary mixture 
O* 'at low density OE]. Recently, the shear viscosity has been also calculated for a dense binary mixture 0. In all these 

cases, the standard procedure to get the transport coefficients is the Chapman-Enskog expansion !& adapted to the 

case of inelastic collisions. 

In the case of elastic collisions, the Chapman-Enskog solution is obtained as an expansion around the local equi- 
^ Librium distribution, while for inelastic collisions the reference state is a local cooling solution with a monotonically 
• <-h jtime decreasing temperature. For spatially homogeneous states, the latter is referred to as the homogeneous cooling 
^^•state (HCS). Since the HCS qualifies as a normal solution, all its time dependence is only through the temperature. 

Nevertheless, in spite of its simplicity, no exact solution of the Boltzmann or Enskog equations describing such a state 
C^H'is known even for one-component systems |10llllj . For a granular binary mixture, the set of coupled Enskog equations 
^ used to describe the system also admits a scaling solution in which the time dependence of the distribution functions 
occurs entirely through the temperature T of the mixture. An important result is that the partial temperatures Ti 
(i = 1,2) of each species (measuring their mean kinetic energies) are different, although their cooling rates are equal 
EH • This effect is generic for multicomponent systems and is a consequence of the inelasticity and the mechanical 
'differences of the particles. This result contrasts with previous studies in granular mixtures |14l 1151 ITffl 117) . where the 
equality of the partial temperatures was assumed. The velocity distribution functions for each species are approxi- 
mately determined by using a Sonine polynomial expansion around Maxwellians, which are defined in terms of the 
temperature for that species. The results show that, in general, the non-Maxwellian corrections have small effects 
on the cooling rates and on the temperature ratio even for strong dissipation. However, the corresponding reference 
Maxwellians for the two species are quite different due to temperature differences. 

As explained above, one of the main difficulties in obtaining the transport coefficients lies in the fact that, in contrast 
to what happens for elastic fluids, the reference state (HCS) depends on time due to the dissipation of energy through 
collisions. To overcome such difficulty, one possibility is to introduce external forces (thermostats) to accelerate the 
particles and hence compensate for collisional cooling. As a consequence, the corresponding reference state is stationary. 
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This mechanism of energy input (different from those in shear flows or flows through vertical pipes) has been used by 
many authors in the past years to analyze different problems, such as non-Gaussian properties (cumulants, high-energy 
tails) of the velocity distribution function long-range correlations an d collisional statistics and short-scale 

structure Since the latter requires the solution of the corresponding linearized hydrodynamic equations around 
the homogeneous state, the explicit expressions for the transport coefficients are needed. The use of Gauss's principle 
of least constraint leads to the Gaussian thermostat, a force proportional to the particle velocity. For the sake of 
brevity, the resulting stationary state will be called here heated stationary state (HSS). It must be pointed out that 
the Enskog equations in the HSS and in the HCS are formally identical when one scales the velocity with the thermal 
velocity [111 1121 IT3] . Therefore, the results obtained for the distribution functions and the temperature ratio 7 in 
the HCS applies to this thermostatted case as well. However, the response of the granular fluid when small spatial 
gradients are introduced is different in both cases. In particular, the shear viscosity calculated in the HCS does not 
coincide with that obtained in the HSS |HJ|H1- The Navier-Stokes transport coefficients of a granular monocomponent 
gas have been obtained by solving the Enskog equation from the Chapman-Enskog method around the HSS 0. Also, 
the shear viscosity of a binary mixture in the HSS has been calculated for dilute 2] and moderately dense [B] systems. 

The partial temperatures and the velocity distribution functions are approximately determined in both the HCS 
and the HSS from a Sonine polynomial expansion Once the HCS and the HSS are well characterized, the 

transport coefficients are given in terms of the solution of linear integral equations. These equations are solved by 
considering the leading terms in a Sonine polynomial expansion as well. A natural question is whether this type of 
approximation leads to accurate results. The direct simulation Monte Carlo (DSMC) method [201 has demonstrated 
to be an efficient and accurate procedure to analyse molecular gases far from equilibrium. This method, conveniently 
adapted to deal with inelastic collisions and finite density systems, provides an alternative route to study the behaviour 
of granular mixtures on the basis of the Boltzmann and Enskog equations, and so to establish the range of validity 
of the approximations used in the theoretical predictions mentioned above. The shear viscosity of a low-density 
granular gas in the HCS was determined from the DSMC method |21|. The results showed a very good agreement 
with the predictions based on the Boltzmann equation in the first Sonine approximation. The numerical experiment 
consisted in preparing an initial in-homogeneous nonequilibrium state corresponding to a transverse shear wave, and 
then analyzing its subsequent evolution in time. The shear wave decayed exponentially with a time scale inversely 
proportional to the viscosity. An alternative route for measuring the shear viscosity consists of preparing a state of 
uniform shear flow (USF) using Lees-Edwards boundary conditions |22|. Macroscopically, this state is characterized by 
constant density, a uniform temperature, and a linear velocity profile. In a ordinary fluid, unless a thermostatting force 
is introduced, the temperature increases in time due to viscous heating. The corresponding energy balance equation 
can be used to determine the shear viscosity for sufficiently long times [221 12H • In a granular fluid, the relationship 
between the temperature and the shear viscosity is not so simple since there is a competition between viscous heating 
and collisional cooling. However, if external forces of the Gaussian form that exactly compensate for the collisional 
energy loss are introduced, the viscous heating effect is still able to heat the system. After a transient period, the 
system reaches a linear hydrodynamic regime where the shear viscosity can be calculated from the pressure tensor. 
This viscosity characterizes the momentum transport when a small perturbation is introduced in the HSS 0IE1- 

The goal of this contribution is to present an overview of the results obtained by means of the Monte Carlo 
simulation for a binary granular mixture in both the HCS and the HSS. The numerical data for the temperature ratio, 
the velocity distribution functions, and the shear viscosity (in the case of the HSS) are compared with the analytical 
results showing an excellent agreement for the range of parameters considered. The plan of the paper is as follows. 
The theoretical framework and the procedure used to calculate the viscosity are described in Sec. |3 The details of 
the simulation method are given in Sec. [31 The analytical results and the simulation data are compared in Sec. 01 We 
close the paper in Sec. with a short summary and conclusions. 

2 Theoretical framework 

Consider a binary mixture of smooth hard spheres of masses mi and ni2 and diameters o~\ and 02- The inelasticity 
of collisions among all pairs is characterized by three independent constant coefficients of normal restitution an, 022, 
and oi\i — a2i, where ay is the restitution coefficient for collisions between particles of species i and j. The molar 
fractions xi = rii/n (i = 1, 2) indicate the relative concentration of the two species, while the density of the mixture can 
be measured by the volume packing fraction <\> — irn/6 (xiaf + a^erf). Here, rii is the number density corresponding 
to species i and n = ri\ + n<i. Of course, the monocomponent case is recovered when one considers mechanically 
equivalent particles, i.e., m\ — rri2, o~\ — 02 and = a. 
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As mentioned in the Introduction, the viscosity characterizing the momentum transport in the HSS is calculated 
in our simulations by considering the USF with a thermostat force that exactly compensates for collisional cooling. 
The shear viscosity can be obtained from the ratio of the rry-element of the pressure tensor to the shear rate for 
sufficiently long times. This viscosity characterizes the evolution of the system when small spatial gradients of the 
velocity field (for instance, a transverse shear wave) are introduced in the HSS. Let us now introduce the Enskog 
equation under USF. At a microscopic level, the USF is generated by Lees-Edwards boundary conditions [221 m which 
are simply periodic boundary conditions in the local Lagrangian frame V = v — a • r and R = r — a • rt. Here, a is the 
tensor with elements a a p = a8 ax 5p v , and a is the shear rate. Using the above frame of reference, the USF becomes 
homogeneous. In terms of these variables, the velocity distribution functions /j(V, t) (i = 1,2) describing the granular 
fluid obey the equations 

8 2 

dtfi - aVy W fi + Tfr = ]T Jf 3 [V\fi(t), f^t)} . (1) 

x i=i 

Here, T is an operator representing the effect of an external force (if it exists) , and jf- is the Enskog collision operator 
Jfi [Vilfi, fj] = o? jXij J dV 2 J da Q(* • g)(<r • g) [a^MV'^f^t) - / 4 (V l7 i)/,(V 2 , t)} . (2) 

In Eq. (0), Xij is the pair correlation function for particles of type i and j when they are at contact, separated by 
Uij = [g{ + <Tj) /2. Note that we have taken into account that Xij is uniform in all the states considered. In this case, 
the Carnahan-Starling approximation provides accurate results. It is given by 

XlJ i - <j> ^ 2 (l - 4>y 2 (l - 4>f ' 1 1 

where 4> = (7rn/6) aijofj and ay = auOjj / '<Jij- Also, in Eq. J2J) cr is a unit vector along their line of centers, 9 is 
the Heaviside step function, g = Vi — V 2 a ■ cr, and 

VJ = Vi - Ni (1 + a- 1 ) (a • g)c? , (4) 



V 2 = V 2 +(i ij (l + aT j 1 )(a-g)a 

+2a a t] d v x , (5) 

where fiij — mi/ [rrii + nrij ) . The collision operators J2J conserve the particle number of each species and the total 
momentum, but the total energy is not conserved. The cooling rate £ is a measure of the loss of energy of the mixture 
due to collisions. In terms of the velocity distribution functions, it is given by 

2 2 

* = iEE^^4(!-4) / rfV i / dV 2 / daeifr-gXa-sffdVirff^t). (6) 

i—l j—1 1 J 

First, let us consider a system that evolves in the absence of shear. Obviously, the analysis can be also carried 
out from the above equations by simply setting a = 0. In the case of elastic particles {puj = 1) and in the absence of 
external forcing [T — 0), it is well-known that the long-time solution of Eq. Q is the Maxwell-Boltzmann equilibrium 
distribution function. On the other hand, if the particles are inelastic (ay < 1) and T = 0, a steady state is not 
possible in uniform situations since, due to the dissipation of energy through collisions, the temperature T decreases 
monotonically with time. In this case, the solution of the Enskog equation tends to the HCS, characterized by the fact 
that all the time-dependence of fi occurs only through the temperature T(t) of the mixture. Dimensional analysis 
requires that fi must be of the form 

Uiv-^^n^m^v*) , (7) 

where vo(t) = y / '2T(t)(mi + 7712) j l {rn\m-i) is a thermal velocity defined in terms of the temperature of the mixture 
T(t), and v* = v/vo(t). Since the time dependence of fi occurs only through T(t), it follows that the temperature 
ratio 7 = T1/T2 must reach a constant value in the long time limit. 

However, by driving a granular fluid by boundaries or external fields it can reach a steady state. The energy 
injected in the fluid may exactly compensate for the energy dissipated by collisional cooling. The same effect can be 
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obtained by means of external forces (thermostats), F' h , acting locally on each particle. These forces are represented 
by the operator T in Eq. (JTJ. Several types of thermostat can be used. Here we consider the so-called Gaussian 
thermostat, i.e., 

Ff = im.CV, (8) 



where C is a positive constant. In this case, 



rfi = c£ r \yfi(y)]- (9) 



The steady state achieved by introducing this thermostat is the HSS. It must be pointed out that the Enskog equation 
in the HSS is formally identical with the Enskog equation in the HCS when one scales the velocity with the thermal 
velocity vo [111 1121 H^] . Therefore, the results obtained for the distribution functions and the temperature ratio 7 
in the HCS applies to this thermostatted case as well. The same does not occur with the shear viscosity. 

In order to calculate the shear viscosity in the HSS, we assume now that the system is under USF. Let us consider 
first the case of elastic collisions. In the absence of a thermostatting force, the energy balance equation yields the 
heating equation [231 EH 

d t T = ~aP xy , (10) 

where P xy is the xy element of the stress tensor P. Since the temperature T increases in time, so does the collision 
frequency v{t) cx WT(t). As a consequence, the reduced shear rate a*(t) = a/v(t) (which is the relevant uniformity 
parameter) monotonically decreases in time and the system asymptotically tends towards that of equilibrium. This 
implies that for sufficiently long times (which means here a* < 1) the system reaches a regime described by linear 
hydrodynamics, and the Navier-Stokes shear viscosity r\ can be identified as j2E] 

v = _ l im . (11) 

This route has been shown to be quite efficient to measure the Navier-Stokes shear viscosity coefficient in ordinary 
fluids [23 [231 123 For a granular mixture, unless a thermostat is introduced, the energy balance equation leads to 
a steady state when the viscous heating effect is exactly balanced by the collisional cooling P [2B1 [2H1 [23 El G2] • 
However, if the granular mixture is excited by the Gaussian thermostat that exactly compensates for the collisional 
energy loss (in this case £ is the instantaneous value of the cooling rate £), the viscous heating still heats the system 
and Eq. (|10|l remains valid. Consequently, the linear relationship (jll|) allows one to determine the shear viscosity 
coefficient in the long time limit. It must be noted that the value of rj calculated in this way corresponds to the 
Navier-Stokes shear viscosity of an excited granular mixture (in the HSS) 00 El, and thus it does not necessarily 
coincide with the value obtained in the unforced case (in the HCS) EH • 

In summary, the Enskog equations Q provide an adequate theoretical framework to study the granular properties 
in the HCS and in the HSS. If one takes a = and T = 0, then the system evolves to the HCS independently of the 
initial conditions considered. For a = and T given by ©, the HSS is achieved. In order to obtain the shear viscosity 
coefficient 77 of the mixture in this state, an arbitrary value o^O has to be taken (USF) and the thermostat force © 
has to be introduced with £ equal to the instantaneous value of the cooling rate £. The coefficient r\ can be measured 
in the long time limit. 

As mentioned in the Introduction, in this context the theoretical predictions are determined by means of a Sonine 
polynomial expansion. In the next section a Monte Carlo algorithm devised to solve JQ) is described. The numerical 
results obtained are compared with the analytical ones in Sec. 01 



3 Monte Carlo simulation method 

The Enskog Simulation Monte Carlo (ESMC) method [531121] is an extension of the well-known direct simulation Monte 
Carlo (DSMC) method [201 to dense gases, and it was devised to mimic the dynamics involved in the Enskog collision 
term. Here we briefly describe the adaptation of this procedure to numerically solve the problem JQ). It must be pointed 
out that since the sought solutions are spatially homogeneous (in the case of the USF thanks to the Lagrangian frame 
used) , the simulation method becomes especially easy to implement and efficient from a computational point of view. 
This is an important advantage with respect to molecular dynamics simulations. Nevertheless, the restriction to this 
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homogeneous states prevents us from analyzing the possible instability of the USF or the formation of clusters or 
microstructures. 

The ESMC method as applied to (IJ is as follows. The velocity distribution function of the species i is represented 
by the peculiar velocities {Vfc} of Ni "simulated" particles: 



1 Ni 

fi(y,t)^ ni —j2 s( y- v k(t)). (12) 



fc=i 

Note that the number of particles N must be taken according to the relation N1/N2 = n\/ri2- At the initial state, 
one assigns velocities to the particles drawn from the Maxwell-Boltzmann probability distribution: 

/i(V, 0) = m V^iO) exp (-V 2 /V£(0)) , (13) 

where V^(0) = 2T(0)/mj and T(0) is the initial temperature. To enforce a vanishing initial total momentum, the 
velocity of every particle is subsequently subtracted by the amount -ZV" 1 Vfc(0). In the ESMC method, the free 
motion and the collisions are uncoupled over a time step At which is small compared with the mean free time and the 
inverse reduced shear rate a* -1 (if a 7^ 0). If those quantities vary significantly, the value of At must be conveniently 
updated in the course of the simulation. 

If a 0, in the local Lagrangian frame particles of each species (i = 1,2) are subjected to the action of a non- 
conservative inertial force F* = —rrii a ■ V. Thus, the free motion stage consists of making Vfc — > Vfc — a • VfcAi. 
In the collision stage, binary interactions between particles of species i and j must be considered. To simulate the 

collisions between particles of species i with j a sample of ^NiU^J^At pairs is chosen at random with equiprobability. 

Hi) 

Here, Wmi is an upper bound estimate of the probability that a particle of the species i collides with a particle of the 
species j. Let us consider a pair {k,£} belonging to this sample. Hereafter, k denotes a particle of species i and I a 
particle of species j. For each pair {k,£} with velocities {Vfc, V^}, the following steps are taken: (1) a given direction 
<Tki is chosen at random with equiprobability; (2) the collision between particles k and I is accepted with a probability 
equal to 0(g k e • vu)J^ r/wmlx, where d$ = 47rcr?.n,-|g w ■ &ki\ and g ke = V fc - V e - try a ■ au\ (3) if the collision is 
accepted, postcollisional velocities are assigned to both particles according to the counterparts of the rule 

Vfc — > Vfc - fj,ji(l + aij)(gkt ■ (tu)&U , (14) 

Vt —> V £ + + Oij)(gki ■ vm)<*u ■ (15) 

If in a collision uj^f > Wmalc, the estimate of Wmax is updated as Wmal = . The procedure described above is 
performed for i = 1, 2 and j = 1,2. 

In the HSS case, after the collisions have been calculated, the thermostat force ||HJl is considered by making 

Vfc ^ Vfc + 1/2 CVfcAi . (16) 

To obtain the shear viscosity value in this state, an arbitrary value o^O has to be taken (USF), and the constant £ in 
lfTH|l must be equal to the instantaneous value of the cooling rate £. The latter is obtained by calculating the granular 
temperature before and after the collision stage. To get the rest of quantities in the HSS, a is set equal to zero and £ 
can be chosen arbitrarily. In fact, this choice only affects to the value of the granular temperature in the steady state. 

In the course of the simulations, one evaluates the distribution function of the species i, its moments and the 
partial temperatures Tj in the usual way. In addition, the kinetic and collisional transfer contributions to the pressure 
tensor are calculated by the expressions 

2 Ni 



^ N 

i=l ' k=l 



]T V fc V *< ( 17 ) 



t 

n \ - 

2NAt 



(j,ijmj(Tij(l + Oij)(gk£ ■ Vki)vkiVki , (18) 



where the dagger means that the summation is restricted to the accepted collisions. The shear viscosity r\ is obtained 
from Ijllfl . while its kinetic contribution is if~ = — lim^oo P^ y /a. In our simulations we have typically taken a total 
number of particles N = Ni + N 2 = 10 5 and a time step At — 3 x lCT 3 Aii/Vbi(0), where An = (V^miOn) -1 is 
the mean free path for collisions 1-1. To improve the statistics, the results are averaged over a number Af = 10 of 
independent realizations or replicas. 
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4 Results 



In this section we compare the analytical predictions mentioned in the Introduction with the numerical results obtained 
from the simulation procedure described above. The dimensionless quantities calculated in both theory and simulation 
are functions of the set of parameters {/i = roi/ro2, u> = <Ji/(?2 1 S = ni/ri2, ctij, <j)}. For the sake of simplicity, in what 
follows we shall assume that an = a 2 2 = ai2 = ot. 

Let us consider first the analysis in the absence of shear rate. The basic quantity measuring the deviation of 
the distribution functions from the corresponding Maxwellians are the cumulants Ci In Fig. 2] we show the 

dependence of c\ and C2 on a for /j, = 2, u) = l, 5=1 and = 0. We also present the corresponding result for 
the one-component system (mechanically equivalent particles, i.e., fi = u> = 1). Note that the results apply in both 
the HCS and in the HSS. The agreement between the simulation data and the theoretical predictions is excellent. 
The small values of these coefficients support the assumption of a low-order truncation in polynomial expansion and 
indicate that the distribution functions for thermal velocities are well represented by the Sonine approximation. To 
confirm this, we have measured the deviation of the distribution functions $>i given by J7J from the corresponding 
Maxwellian. More precisely, we have evaluated the function Ai(v*) defined by the relation 



*<(«*) 



3/2 



-XiV* 



(19) 



where A, ; = Tj (Tj/x^). The function Ai(v*) is plotted in Fig.Hfor n = 4, u> = 1, 6 = 1/2, a = 0.5 and = 0. The 
dashed line is the first Sonine approximation 
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(20) 



As could be expected, the simulation curve agrees very well with the corresponding Sonine polynomial, confirming the 
accuracy of the analytical solution in the region of thermal velocities. 



0.15 



0.10 



0.05 



0.00 



-0.05 




_l i L 



0.5 0.6 0.7 0.8 0.9 1.0 

a 



Figure 1: Plot of the coefficients Cj versus the restitution coefficient a for fj, = 2, u> = 6 = 1 and = 0. The solid line 
and the circles refer to ci while the dashed line and the squares correspond to c 2 . The dotted line and the triangles 
refer to the common value in the single component case. The lines are the theoretical predictions and the symbols 
correspond to the simulation results. 

One of the main new results of the description showed in Refs. ^JEj is that the partial temperatures are different 
(7 7^ 1). This conclusion contrasts with previous results derived for granular mixtures |14l 1151 IT^l I17j . where it was 
implicitly assumed the equipartition of granular energy between both species (i.e., 7 = 1). In Figs. [2 ^ and [S] we 
plot the temperature ratio 7 versus the restitution coefficient a for different choices of the mechanical parameters 
characterizing the mixture. The theory as well as the simulation results clearly indicate that 7 is different from unity, 
even for weak inelasticity. Figure shows the dependence of 7 on a for w = l,<5 = 2,0 = O and several values of 
the mass ratio. The agreement between theory and simulation is very good, implying the accuracy of the expression 
of 7 obtained in the first Sonine approximation. For large differences in the mass ratio, the temperature differences 
are significant. As can be observed in Fig. the influence of the concentration ratio 6 on the temperature ratio 7 is 
not as strong as that observed with the mass ratio, although is still quite important. The dependence of the relative 
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Figure 2: Plot of the simulation values of the function A\(v*) defined by Eq. I|19f> for fi = 4, u> = 1, 6 = 1/2, a = 0.5 
and (f> = 0. The dashed line is the Sonine polynomial J3UJ. 

temperature ratio 7(0?, (j))/j(a, 0) on the volume packing fraction (f> is plotted in Fig. for /j = 2, w = 2, 6 = 1/2, 
and two different values of a: a = 0.6 and 0.8. We see that, for a given value of the density, the relative temperature 
ratio decreases as the degree of inelasticity increases. It is evident again the excellent agreement between the Sonine 
predictions and the simulation data. 
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Figure 3: Plot of the temperature ratio 7 versus the restitution coefficient a for (j = l ) 8 = 2 ) <f> = and three different 
values of the mass ratio: /.i = 1/10 (dotted line and triangles), fi ~ 2 (dashed line and squares) and /i = 10 (solid line 
and circles). The lines are the theoretical predictions and the symbols correspond to the simulation results. 

Let us focus on the sudy of the shear viscosity coefficient r\. The theoretical analysis carried out in Refs. [3 IE] 
shows that the transport properties are affected by the thermostat introduced. In particular, the values obtained 
for the shear viscosity in the HCS and in the HSS are different. In Fig. we plot the reduced shear viscosity 77/770, 
770 being the corresponding elastic value, in the HCS and in the HSS for a monocomponent gas. In both cases the 
viscosity increases with the dissipation. However, at a quantitative level, the influence of dissipation on the viscosity 
in the HCS case is smaller than in the HSS. The agreement between theory and simulation is quite good. The same 
trends are observed for a binary mixture (cf. Fig. 1 in Ref. 8 ). 

Next, in Figs.[7||Hland|Hlwe show the influence of dissipation on the reduced shear viscosity rj(a)/r]o for a heated 
granular binary mixture at low density ((f) — 0) with different values of the mass ratio, the size ratio, and the 
concentration ratio. Three different values of the restitution coefficient are considered: a = 0.9,0.8, and 0.7. In Fig. H 
we plot the ratio r/(a)/r]o versus the mass ratio /x for u> = 8 = 1. Again, the symbols represent the simulation data while 
the lines refer to the theoretical results obtained from the Boltzmann equation in the first Sonine approximation. We 
see that in general the deviation of 77(a) from its functional form for elastic collisions is quite important. This tendency 
becomes more significant as the mass disparity increases. The agreement between the first Sonine approximation and 
simulation is seen to be in general quite good. This agreement is similar to the one found in the monocomponent case 
(Fig.EJ. At a quantitative level, the discrepancies between theory and simulation tend to increase as the restitution 
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Figure 4: Plot of the temperature ratio versus the restitution coefficient a for fx = 4, u> = 1, <f> = and three different 
values of the concentration ratio: 6 = 1/4 (dotted line and triangles), 5 — 1 (dashed line and squares) and 5 = 4 (solid 
line and circles). The lines are the theoretical predictions and the symbols correspond to the simulation results. 



Figure 5: Plot of the relative temperature ratio 7(a, (j>)/~f(a, 0) as a function of the volume packing fraction <f> for 
fi = oj = 2, 8 = 1/2, and two different values of the restitution coefficient: a = 0.6 (solid line and circles) and a = 0.8 
(dashed line and squares). The lines are the theoretical predictions and the symbols correspond to the simulation 
results. 

coefficient decreases, although these differences are quite small (say, for instance, around 2% at a — 0.7 in the disparate 
mass case toi/to2 = 10). 

The influence of the size ratio on the shear viscosity is shown in Fig. [S] for /i = 4 and 6—1. We observe again 
a strong dependence of the shear viscosity on dissipation. However, for a given value of a, the influence of lu on r\ 
is weaker than the one found before in Fig. 0for the mass ratio. The agreement for both a = 0.9 and a = 0.8 is 
quite good, except for the largest size ratio at a = 0.8. These discrepancies become more significant as the dissipation 
increases (say, for instance, a = 0.7), especially for mixtures of particles of very different sizes. Finally, Fig. shows 
the dependence of rj(a)/rjo on the concentration ratio for fi = 4 and u> = 1. We observe that both the theory and 
simulation predict a very weak influence of composition on the shear viscosity. With respect to the influence of 
dissipation, the trends are similar to those found in Figs. and |H1 the main effect of inelasticity in collisions is to 
enhance the momentum transport with respect to the case of elastic collisions. The agreement now between theory 
and simulation is very good, even for disparate values of the concentration ratio and/or strong dissipation. Therefore, 
according to the comparison carried out in Figs. El [HI and|§| we can conclude that the agreement extends over a wide 
range of values of the restitution coefficient, indicating the reliability of the first Sonine approximation for describing 
granular flows beyond the quasielastic limit. 

In order to analyze density effects on the shear viscosity, in Figs . 1 1 01 and 1 1 II we plot the kinetic part and the total 
value of this quantity versus the volume packing fraction. The parameters of the mixture are fj, = 4, u> = 1, and 
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Figure 6: Reduced shear viscosity 77/ 770 as a function of the restitution coefficient a for a monocomponent gas at low 
density (0 = 0). The solid lines refer to the analytical results derived in the HCS (a) and in the HSS (b). The symbols 
are the results obtained from the simulation in the HSS. 
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Figure 7: Plot of the ratio 77/770 as a function of the mass ratio /i for uj = S = 1 and three different values of the 
restitution coefficient a: (a) a = 0.9 (circles), (b) a — 0.8 (squares), and (c) a — 0.7 (triangles). The lines are the 
theoretical predictions and the symbols refer to the simulations. 

6=1. Three different values of a are studied: a = 0.9, 0.8, and 0.7. Figure [Till shows the dependence of the kinetic 
part r) k * = vrj k /nT on the solid fraction 0, while the total shear viscosity 77* = vr\l nT is plotted in Fig. ^] Here, 
v = s/nna^VQ is an effective collision frequency. The good agreement between theory and simulation indicates that 
both kinetic and collisional transfer contributions are given accurately by the first Sonine approximation. As in the 
monocomponent case (cf. Fig. 2 in Ref. [Hj), the shear viscosity of a granular mixture decreases (increases) as the 
inelasticity increases if the solid fraction is larger (smaller) than a given threshold value 0o. The value of 0o depends 
on the parameters of the mixture although it is practically independent of dissipation. For the mixture considered in 
Fig.EJJ 0o (a) ^ 0.22. 

Finally, we explore in Fig. 1121 the influence of dissipation on the reduced shear viscosity 77* for different values 
of the mass ratio in a mixture at finite density (0 = 0.2). As in the low-density case, we see that the influence of 
dissipation on 77* becomes important as the mass disparity increases. At a given value of the mass ratio, 77* decreases 
(increases) with dissipation if the mass ratio is smaller (larger) than a certain threshold value, the value of which 
seems to be again practically independent of the restitution coefficient. Regarding the comparison between kinetic 
theory and simulation, we see that the agreement between both approaches is similar to the one previously obtained, 
although the discrepancies tend to increase as a decreases. 
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Figure 8: Plot of the ratio 77/770 as a function of the size ratio oj for // = 4, 6 = 1 and three different values of the 
restitution coefficient a: (a) a = 0.9 (circles), (b) a = 0.8 (squares), and (c) a — 0.7 (triangles). The lines are the 
theoretical predictions and the symbols refer to the simulations. 
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Figure 9: Plot of the ratio 77/770 as a function of the concentration ratio S for fj, = 4, u = 1 and three different values 
of the restitution coefficient a: (a) a = 0.9 (circles), (b) a = 0.8 (squares), and (c) a — 0.7 (triangles). The lines are 
the theoretical predictions and the symbols refer to the simulations. 

5 Concluding remarks 

In this contribution we have presented an overview of recent results on Monte Carlo simulations of a granular binary 
mixture. We have considered two homogeneous states: the so-called homogeneous cooling state (HCS) and the 
stationary state achieved when a gaussian thermostat force is introduced (HSS). It is well known that the results 
obtained for the velocity distribution functions in both cases are identical when they and velocities are conveniently 
scaled [njEIIEl- Transport properties (such as the shear viscosity coefficient) are related to the response of the 
granular fluid when small spatial gradients are introduced. In this case, the external forces used as thermostats do 
not play a neutral role. For instance, Fig. shows that the discrepancy between the shear viscosity of the driven and 
undriven cases is quite significant. It must be noticed that the above conclusion only affects to this type of external 
forcing mechanisms, since driving the system by shaking, vibration, and even the action of a weak external field (such 
as the gravity field) does not modify the transport coefficients of the granular fluid. 

The main motivation of performing Monte Carlo simulations was to verify the range of validity of the approximate 
analytical predictions obtained from the Enskog equation by using a Sonine polynomial expansion in both the HCS 
and the HSS. In general, the agreement found is excellent in all the cases considered, even for large disparities of 
masses, sizes and concentrations. Discrepancies slightly increase as energy dissipation in collisions increases. This 
agreement is a further testimony to the validity of a hydrodynamic description for granular media beyond the weak 
dissipation limit. 





10 



1.3 
1.2 
1.1 

k» 

V 1.0 
0.9 
0.8 
0.7 
0.6 

0.000 0.075 0.150 0.225 0.300 

Figure 10: Plot of the kinetic part rj k * of the reduced shear viscosity as a function of the solid fraction for /j, = 4, 
to = S = 1 and three different values of the restitution coefficient a: a — 0.9 (solid line and circles), a = 0.8 (dashed 
line and squares), and a = 0.7 (dotted line and triangles). The lines are the theoretical predictions and the symbols 
refer to the simulations. 
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Figure 11: Plot of the reduced shear viscosity ry* as a function of the solid fraction <fi for /i = 4, lj = 5 = 1 and three 
different values of the restitution coefficient a: a — 0.9 (solid line and circles), a = 0.8 (dashed line and squares), and 
a = 0.7 (dotted line and triangles). The lines are the theoretical predictions and the symbols refer to the simulations. 

Nevertheless, the analysis based on the Enskog equation (JJJ is limited in several aspects. A test of the utility of the 
Enskog equation at high densities is possible using molecular dynamics simulations. Previous comparisons at the level 
of partial temperatures jSHJ and self-diffusion coefficient [121 indicate that the range of densities for which the Enskog 
theory applies decreases with increasing dissipation. It must be noted that upon describing the HCS and the HSS by the 
Enskog equations , it has been implicitly assumed the validity of the "molecular chaos" hypothesis of uncorrelated 
binary collisions. However, molecular dynamics simulations of hard disks have shown a non-uniform distribution of 
impact parameters for high enough dissipation |.'i5| . In addition, there exist long range spatial correlations in density 
and flow fields which can not be understood on the basis of the Enskog equation ISS] • These two effects are associated 
with the appearance of the so-called cluster instability [67\ for systems sufficiently large. Since we have simulated 
directly the spatially uniform equation such an instability is precluded in the simulations. 

The main finding of the analysis carried out for binary mixtures is that, in general, the partial temperatures of 
each species are different ^2 El, in contrast to what was assumed in previous studies [141 1151 H?>l 117) . This violation 
of energy equipartition has been subsequently confirmed in experiments of vibrated granular mixtures [HBI EH! an d m 
recent molecular dynamics simulations of the HCS [HSj- The Monte Carlo simulations collected in this contribution, 
and other ones not shown here, demonstrate that the analytical predictions accurately capture the dependence of the 
temperature ratio on the parameters of the mixture. 

The violation of energy equipartition has been also found in granular mixtures under simple shear flow |29l 1301 1311 
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Figure 12: Plot of the reduced shear viscosity -q* as a function of the mass ratio /1, for u = S = 1, <j> = 0.2 and three 
different values of the restitution coefficient a: a = 0.9 (solid line and circles), a = 0.8 (dashed line and squares), and 
a = 0.7 (dotted line and triangles). The lines are the theoretical predictions and the symbols refer to the simulations. 

132) . The simple shear flow state has been extensively studied for molecular fluids as a prototype problem to analyse 
nonlinear transport properties. Nevertheless, the nature of this state is quite different in the case of granular fluids. 
While for clastic fluids the temperature grows monotonically in time due to viscous heating [531 [^J, a steady state is 
possible for granular media when the effect of viscosity is exactly compensated by the collisional cooling. In this case, 
the system reaches a steady state and the temperature achieves a constant value. This state can be analysed using 
the algorithm described in Sec.|3|by taking an arbitrary value of the shear rate a and setting £ = 0. The value of the 
reduced shear rate a* — a/v, which is the uniformity parameter that characterizes the steady state, is a function of the 
normal restitution coefficients ay. The results presented in Refs. |29l 1301 l3"Tl 132) show again that the temperature ratio 
is clearly different from unity (as may be expected since the system is out of equilibrium) and strongly depends on the 
restitution coefficients as well as on the parameters of the mixture. The influence of the temperature differences on 
the rheological properties was analyzed from the Boltzmann kinetic theory by using a Sonine polynomial expansion. 
The comparison between kinetic theory and Monte Carlo simulations showed an excellent agreement over the range 
of parameters investigated. 

Partial support from the MCYT (Spain) through Grant No. ESP2003-02859 is acknowledged. 
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